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I. Summary 

A class of explicit and implicit total variation diminishing (TVD) schemes for the compressible 
Euler and Navier-Stokcs equations has been developed 1 1-4). They have the property of not 
generating spurious oscillations across shocks and contact discontinuities. In general, shocks 
can be captured within 1-2 grid points. For the inviscid case, one can divide these schemes into 
upwind TVD schemes and symmetric (non-upwind) TVD schemes. The upwind TVD scheme is 
based on the second-order TVD scheme developed by Ilarten (5). The symmetric TVD scheme 
developed by the author is a generalization of Roe [6] and Davis’s [7] TVD Lax-WendrofT scheme. 
The objective of this paper is to investigate the performance of these schemes on some viscous 
and inviscid airfoil steady-state calculations. A comparison of the symmetric and upwind TVD 
schemes is included. 

II. Description of Algorithm for System of Hyperbolic Conservation Laws 

The notion of upwind and symmetric TVD schemes, including formulation and extension to 
system cases (in uniform Cartesian grids), can be found in references (1-3,5) . Here the extension 
of the implicit second-order-accurate TVD scheme for hyperbolic systems of conservation laws 
in curvilinear coordinates (4) is briefly described. 

Consider a two-dimensional system of hyperbolic conservation laws 

d.i) 

31 ax dy 

Here Q, F[Q) and G(Q) are column vectors of m components. 

A generalized coordinate transformation of the form £ -• £(i,y) and r ; — t]{x,y) which 
maintains the strong conservation law form of equation (2.1) is given by 


dQ dF{Q) 3G(Q) _ 

at + “af“ “aT " ’ 


( 2 . 2 ) 


where Q - Q/J, F - (£ X F I ( V G)/J, G =• (i j x F + 1 i y G)/J, and J ---• (, T r h - the Jacobian 
transformation. Let A ~ dF/BQ and B ~ BG/3Q\ then the Jacobians A and B of F and (» 
can be written as 

A ~ (£ X A + £ V D) (2.3a) 

D = (r/ z A + rj v B). (2.3b) 

Let the eigenvalues of A be (uj,a|, ...,a^) and the eigenvalues of B be (a^,a*,...,o”). Denote 
li^ and R n as the matrices whose columns are eigenvectors of A and B, and denote R ^ 1 and 
R~ 1 as the inverses of R $ and R,,. 

Let the grid spacing be denoted by A£ and At; such that £ ~ jA£ and tj ~ fcAr/. Let i *, 
denote some symmetric average of Qj,k and Qj + i,k (for example, Q J + k = 0.5* (Q } .f ■i,k + Qj,k)> 
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or the Roc’s average |S| for gas dynamics). Leta' J + j, /?“ 4 \ denote the quantity- a' ( , li(, 

Ji ~ l related to A evaluated at Q J( i k . Similarly, let a' k + i, H k{ i, /? M ‘ , denote the quantities 
a[,, 7i n , li n l related to 71 evaluated at Q } ,k+{- 


Define 
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as the difference of the characteristic variables in the local {-direction, and define 
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(2.4b) 


as the difference of the characteristic variables in the local q -direction. The J }ik is the Jacobian 
transformation evaluated at (jA{, A:Ar?). The averaged Jacobians are used here in order to 
preserve frccstrcain. 

With the above notation, a one-parameter family of TVD schemes can be written as 
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(2.5a) 


where 0 < 6 < 1, A f - At/A{, A'' = At/A??, and At is the time step. This one-parameter 
family of schemes contains implicit as well as explicit schemes. When 9 ~ 0, (2.5) is an explicit 
method; when 6 f 0, it is an implicit scheme. For example, if 9 = 1/2, the time differencing 
is the trapezoidal formula and scheme (2.5) is second-order in time and space. If 9 ~ 1, the 
time differencing is the backward Euler method and scheme (2.5) is first-order in time but 
second-order in space. 


The numerical flux function i k for both the upwind and symmetric. TV D schemes [2,3] 
can be expressed as 


F, 


ji i .*= 




(2.5b) 


Similarly, we can define the numerical flux Gf jk ^.^. 

Upwind TVD Scheme [5] 

For a particular form of the upwind TVD scheme [2,5], the elements of the + i denoted by 

(^ + i )-.* = are 


with 
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g‘ } = minmod(a^_ j , a' + * ) . 


(2.6a) 


(2.6b) 


The minmod function of a list of arguments is equal to the smallest number in absolute value if 
the list of arguments is of the same sign, or is equal to zero if any argument is of opposite sign. 
The function V> * s 
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Here i is a small positive parameter (see reference 5 for a formula oft), and 
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whore a 1 . 4 are elements of (2.4a). 
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Symmetric T VI) Scheme |3j 

For two particular forms of the symmetric TVD scheme, the elements of the <J> j+ j 
by (<A‘ + i) S are 


and 


(^») 5 = - V ' K ^)( 2 o ; 4 , 


)' 


denoted 


(2.7a) 

(2-7b) 


A Conservative Linearized ADI Form For Steady-State Application 

For two-dimensional steady-state applications, the implicit schemes (2.5) can be solved by - 
some appropriate relaxation method other than alternating direction implicit (ADI) and will 
be the direction of future research. The schemes considered in this paper are implemented in 
a conservative noniterative ADI form |2|. For steady-state applications, the numerical solution 
is independent of the time stop. The implicit operator has a regular block tridiagonal structure 
and the resulting block tridiagonal matrix is diagonally dominant. One can modify a standard 
central difference code by simply changing the conventional numerical dissipation term into the 
one designed for the TVD scheme; i.e., the third term of equation (2.5b). - The only difference 
in computation is that the current scheme requires a more elaborate dissipation term for the 
explicit operator; no extra computation is required for the implicit operator. For the Navicr- 
Stokes applications, the convection terms are discretized by a TVD scheme, and the diffusion 
terms are discretized by central difference approximations. 

A conservative linearized ADI form of equation (2.5) used mainly for steady-st.at.c applications 
as described in details in reference [2] can be written as 
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is used and 
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Here A j4Uk , are (2.3) evaluated at (j + 1 ,k) and [j,k 4- 1). The expression diag(z‘) 

denotes a diagonal matrix with diagonal elements z l . All of the inviscid calculations shown in 
this paper use (3.1). 
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l-or steady-state application, a simple algorithm utilizing the TVD scheme for the Navier- 
Stokes equations is to difference the hyperbolic terms the same way as before, and then central 
difference the viscous terms. The final algorithm is the same as equation (3.1) except that 
the spatia! central differencing of the viscous term is added to the right hand side of equation 
(3.1. I he numerical solution shown below illustrates that this algorithm produces a fairly good 
solu Hon : for the case of a HA 102822 airfoil calculation. A treatment for time-accurate calculations 
can be found in reference (2). 

Numerical Results 

V°'T»m f ‘T i,1ViS< ; i( ! ralcula, ’ i<ms * u P wi,,d TVD ■‘■•hcn.08 produce sharper shocks than sym- 
metric VI) schemes |7|. For the current two methods (2.6) and (2.7), this seems to be not the 
case. I he two methods appeared to produce almost identical results for flow field conditions 
ranging from subcnt.cal to transonic and supersonic. Figures 1 and 2 show the comparison of 
equation (2.6) with (2.7b) for two inviscid steady-state airfoil calculations. The advantages of 
symmetric I VI) schemes are that they require less computational efTort and provide a more 
natural way of extending the scheme to two and three-dimensional problems. 

Figures 3 and 4 show an inviscid comparison of the symmetric TVD scheme with the widely 
distributed computer code ARC2I), version 150 (9). The freestream Mach numbers are A/ (v l 
and 18, and the angle of attack is « =■ 7. The pressure coefficient distributions (not shown) are 
Jdent.ca between the two methods and yet the flow field appears very different. The symmetric 
1 VD scheme gives a very well-ordered flow structure and can still capture the shocks with a 
coarse grid, especially near the trailing edge of the airfoil. On the other hand, the ARC2D 
code did rather poorly. The same problem was studied for the upwind TVD scheme and the 
results and convergence rates were found to be almost identical to those for the symmetric TVD 
scheme. A resid uai of 10“ can be reached at around 400-506 steps, ARC2D, however, required 
only 200-300 steps to converge to the same residual. The ARC2D, version 150 computer code is 
based on the Beam and Warming ADI algorithm |10| but uses a mixture of second and fourth- 
order numerical dissipation terms. Figure 5 is an example of the viscous case for the RAK2822 
air o. using the upwind TVD scheme. The thin layer Navier-Stokes equations witl, the algebraic 
turbulence model of Baldwin and Lomax |12| are used and the transition is fixed at 3% chord 
1 he overall agreement with experiments is quite good. The i 2 - norm residual of 10" 7 can be 
reached in around 900 steps. 

Concluding Remarks 

Both the symmetric and upwind TVD schemes are designed to capture shock waves accurately 
while not exlnb. ting the spurious oscillation associated with the more classical second-order 
schemes. Numerical experiments with the TVD schemes on problems containing no shock llll 
show that there is no advantage of the TVD scheme over the conventional bax-WcndrofT type 
scheme. Numerical experiments with various inviscid airfoil calculations show that for weak to 
moderate shocks the main advantage of the TVD schemes is that one can capture the shocks 
in 1-2 grid points |5). '1 he flow field away from the shock looks very much like the classical 
second-order central difference methods. But as one increases the shock strength, especially in 
supersonic flow, '] VI) schemes provide superior flow field solutions even with a very coarse yet 
non-clustering grid near shock waves. Numerical experiments also indicate that the symmetric 
and upwind VD schemes are also applicable for viscous calculations. The current study further 
shows that the symmetric TVD scheme is just as accurate as the upwind TVD scheme while 
requiring less computational effort than the upwind TVD scheme. 
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Fig. 3 Comparison of asymmetric TVD (SYMTVD) scheme with ARC2D (version 150) for the Mach contours, 
pressure contours and entropy contours of the NACA0012 airfoil with M oo = 1.2, a - 7 using a 163 x 49 
C grid as shown in figure 1. 
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Fig. 4 Comparison of a symmetric TVD (SYMTVD) scheme with ARC2D (version 150) for the Mach contours 
pressure contours and entropy contours of the NACA0012 airfoil with = 1.8 a = 7 usinrr a 163 x 49 
C gnd as shown in figure 1. * 
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